信号与系统 Mathematica 代码

第 1 章 信号变换与运算

1.1 常用信号

​x
1
(* 
2
   mathematica 代码比 matlab 好写, 而且软件还挺轻便
3
   1 指数
4
   1.1 指数信号
5
*)
6
Plot[E^t, {t, -1, 3}]
7
(* E^t 也可改为 Exp[t] *)
8

9
(* 1.2 单边指数信号 *)
10
Plot[E^t UnitStep[t], {t, -1, 3}]
11
(*
12
    UnitStep[t] 是单位阶跃函数;
13
    单位阶跃函数的输入可以是向量;
14
*)
15

16
Plot[E^t Boole[t > 0], {t, -1, 3}]
17
(* 
18
    简单的表达式用 UnitStep 很方便；
19
    复杂的表达式用 Boole 很方便.
20
*)
21

22
(* 1.3 复指数信号 *)
23
Plot[Re[E^((-1 + 10 I) t)], {t, -1, 3}]
24
Plot[Im[E^((-1 + 10 I) t)], {t, -1, 3}]
25
Plot[Abs[E^((-1 + 10 I) t)], {t, -1, 3}]
26
Plot[Arg[E^((-1 + 10 I) t)], {t, -1, 3}]
27

28
(* 2 正弦 *)
29
(* 2.1 正弦信号 *)
30
Plot[Sin[2*Pi*t + Pi/6], {t, 0, 2 Pi}]
31

32
(* 4.2 单位阶跃信号 *)
33
Plot[
34
    UnitStep[t], {t, -1, 3},
35
    GridLines -> {Automatic, Automatic}
36
] (* 绘制网格线 *)
37
(* 
38
    HeavisideTheta[t] 也可以表示单位阶跃函数
39
    HeavisideLambda[t] 是一个周期的三角分布
40
    HeavisidePi[t] 表示矩形脉冲
41
*)
42

43
(*
44
    4.3 单位冲激信号
45
    理论上是这样的, 不过从图像上看不出来.
46
    也可以用 UnitStep'[t] 表示
47
*)
48
Plot[DiracDelta[t], {t, -1, 3}]

1.2 信号变换

xxxxxxxxxx
5
1
f1[t_] := t^2
2
(* f1[t_] := . 可以取消定义*)
3
f2[t_] := f1[3 - 2 t]
4

5
Plot[f2[t], {t, -1, 3}]

1.3 信号分解

xxxxxxxxxx
11
1
f[t_] := Cos[t] + Sin[t] + 1
2

3
(* 分解为直流分量与交流分量 *)
4
fd[t_] := Integrate[t1 f[t1], {t1, -Pi, Pi}]
5
fa[t_] := f[t] - fd[t]
6
Plot[fd[t], {t, -Pi, Pi}]
7
Plot[fa[t], {t, -Pi, Pi}]
8

9
(* 分解为奇分量与偶分量 *)
10
Plot[(f[t] - f[-t])/2, {t, -Pi, Pi}]
11
Plot[(f[t] + f[-t])/2, {t, -Pi, Pi}]

第 2 章 连续时域分析

2.1 系统响应

xxxxxxxxxx
7
1
(* 似乎无法求解冲激响应, 如直接使用 DiracDelta, 则结果恒为零 *)
2
solution = DSolve[{
3
    2 y''[t] + y'[t] + 8 y[t] == UnitStep[t],
4
    y[0] == 0, y'[0] == 0
5
}, y[t], t]
6
r[t_] := y[t] /. solution[[1]]
7
Plot[r[t], {t, 0, 25}, PlotRange -> All]

2.2 函数卷积

xxxxxxxxxx
9
1
(* 以一个周期的三角波为例, 首先定义信号, 也可以直接使用 HeavisideLambda[t] *)
2
f[t_] := Which[-1 <= t <= 0, 1 + t, 0 < t <= 1, 1 - t]
3

4
(* 法一: 利用卷积定义 *)
5
Integrate[f[x]*f[t - x], {x, -1, 1}]
6
Plot[%, {t, -2, 2}]
7

8
(* 法二: 利用卷积函数 *)
9
Plot[Convolve[f[x], f[x], x, t], {t, -2, 2}]

2.3 卷积的逆

xxxxxxxxxx
18
1
(* 离散的傅里叶变换及其逆变换 *)
2
x = {1, 1, 2, 2, 1, 1, 0, 0}
3
X = Fourier[x]
4
InverseFourier[X] (* 可能有误差 *)
5

6
(* 连续的傅里叶变换及其逆变换 *)
7
F = FourierTransform[Sin[t]/t, t, \[Omega]]
8
f = InverseFourierTransform[F, \[Omega], t]
9

10
(* 已知 f2 和 f 根据 f = f1 * f2 求 f1 *)
11
f1 = UnitStep[t];
12
f2 = t UnitStep[t];
13
f = Convolve[f1, f2, t, x]
14
f1Test = InverseFourierTransform[
15
    FourierTransform[f, x, \[Omega]] /
16
    FourierTransform[f2, x, \[Omega]],
17
    \[Omega], x
18
] (* 理论上应该是这样的, 但是实测无法化简 *)

2.4 可视化卷积

xxxxxxxxxx
23
1
DynamicModule[{f, g, h},
2
 f[x_] := UnitBox[x/2];
3
 g[x_] := UnitBox[x];
4
 h[y_] = Convolve[f[x], g[x], x, y];
5
 Manipulate[
6
  GraphicsRow[{
7
    Plot[{f[x], g[y - x]}, {x, -5/2, 5/2}, Sequence[
8
     Filling -> Axis, PlotRange -> {-1, 1.5}, Exclusions -> None, 
9
      Ticks -> {Automatic, {0, 1}}, AxesLabel -> Automatic, 
10
      PlotLabel -> Style["f(x), g(y-x)", Italic], Epilog -> {
11
Arrowheads[Small], 
12
Arrow[{{y, -0.75}, {y, -0.1}}], 
13
Text[
14
Style["y", Italic], {y + 0.4, -0.75}]}]],
15
    Plot[f[x] g[y - x], {x, -5/2, 5/2}, Sequence[
16
     Filling -> Axis, PlotRange -> {-1, 1.5}, Exclusions -> None, 
17
      Ticks -> {Automatic, {0, 1}}, AxesLabel -> Automatic, 
18
      PlotLabel -> Style["f(x)\[CenterDot]g(y-x)", Italic]]],
19
    Plot[h[x], {x, -5/2, y}, Sequence[
20
     Filling -> Axis, PlotRange -> {{-(5/2), 5/2}, {-1, 1.5}}, 
21
      Ticks -> {Automatic, {0, 1}}, AxesLabel -> {
22
Style["y", Italic], None}, PlotLabel -> Style["(f*g)(y)", Italic]]]}, 
23
   ImageSize -> 480], {{y, -0.8, Style["y", Italic]}, -2, 2}]]

第 3 章 连续频域分析

3.1 周期信号的频谱

3.1.1 周期矩形脉冲

xxxxxxxxxx
20
1
(* 例 3.13: 周期矩形脉冲的双边频谱 (由于结果为实数, 直接绘出频谱) *)
2
A = 1; T = 4;
3
FnPoints = Table[
4
    {n, Integrate[
5
        A HeavisidePi[t] E^(-I*2 n Pi*t/T),
6
        {t, -1, 1}
7
    ] / T},
8
    {n, -20, 20}
9
] (* 不能分开写 Fn = Integrate[...], 否则结果无法代入 n = 0. *)
10
ListPlot[FnPoints, PlotRange -> All, Filling -> Axis]
11

12
(* 或者这样写 *)
13
A = 1; T = 4;
14
FnPoints = Table[
15
    {n, FourierCoefficient[
16
        A HeavisidePi[t], t, n,
17
        FourierParameters -> {1, 2 Pi/T}
18
    ]}, {n, -20, 20}
19
];
20
ListPlot[FnPoints, PlotRange -> All, Filling -> Axis]

3.1.2 周期三角波信号

xxxxxxxxxx
22
1
(* 例 3.14: 周期三角波信号的双边频谱 *)
2
f[t_] := Which[
3
    -1 <= t <= -0.5, -2 t - 2,
4
    -0.5 < t <= 0.5, 2 t,
5
    0.5 < t <= 1, 2 - 2 t
6
] (* 周期三角波 *)
7

8
AbsFnPoints = Table[
9
    {n, Abs[FourierCoefficient[
10
        f[t], t, n,
11
        FourierParameters -> {1, 2 Pi/2}
12
    ]]}, {n, -10, 10}
13
]; (* 傅里叶系数的模长 *)
14
ListPlot[AbsFnPoints, PlotRange -> All, Filling -> Axis]
15

16
ArgFnPoints = Table[
17
    {n, Arg[FourierCoefficient[
18
        f[t], t, n,
19
        FourierParameters -> {1, 2 Pi/2}
20
    ]]}, {n, -10, 10}
21
]; (* 傅里叶系数的幅角 *)
22
ListPlot[ArgFnPoints, PlotRange -> All, Filling -> Axis]

3.2 非周期信号的频谱

3.2.1 单边指数信号

xxxxxxxxxx
7
1
(* 例 3.15: 单边指数信号的傅里叶变换 *)
2
F = FourierTransform[
3
    E^(-3 t) UnitStep[t], t, \[Omega],
4
    FourierParameters -> {1, -1}
5
] (* 傅里叶变换 *)
6
Plot[Abs[F], {\[Omega], -6, 6}, GridLines -> Automatic]
7
Plot[Arg[F], {\[Omega], -6, 6}, GridLines -> Automatic]

3.2.2 门信号的频谱

xxxxxxxxxx
13
1
F = FourierTransform[
2
    HeavisidePi[t/2], t, \[Omega],
3
    FourierParameters -> {1, -1}
4
]
5
Plot[
6
    HeavisidePi[t/2], {t, -2, 2},
7
    GridLines -> Automatic
8
] (* 门信号 *)
9
Plot[
10
    F, {\[Omega], -30, 30},
11
    GridLines -> Automatic,
12
    PlotRange -> All
13
] (* 门信号的频谱 *)

3.3 傅里叶变换的性质

3.3.1 对称性质

xxxxxxxxxx
15
1
Plot[
2
    FourierTransform[
3
        Sinc[t]/Pi, t, \[Omega],
4
        FourierParameters -> {1, -1}
5
    ], {\[Omega], -2, 2},
6
    GridLines -> Automatic
7
] (* Sa 函数的傅里叶变换 *)
8
Plot[
9
    FourierTransform[
10
        HeavisidePi[t/2]/2, t, \[Omega], 
11
        FourierParameters -> {1, -1}
12
    ], {\[Omega], -20, 20}, 
13
    GridLines -> Automatic,
14
    PlotRange -> All
15
] (* 门函数的傅里叶变换 *)

3.3.2 卷积定理

xxxxxxxxxx
25
1
Plot[
2
    FourierTransform[
3
        HeavisideLambda[2 t], t, \[Omega], 
4
        FourierParameters -> {1, -1}
5
    ], {\[Omega], -30, 30}, 
6
    GridLines -> Automatic,
7
    PlotRange -> All
8
] (* 直接计算三角波信号的傅里叶变换 *)
9

10
Plot[
11
    Convolve[2HeavisidePi[2x], HeavisidePi[2x], x, t],
12
    {t, -2, 2}
13
] (* 验证门函数的卷积是三角波 *)
14

15
Plot[
16
    2 FourierTransform[
17
        HeavisidePi[2t], t, \[Omega], 
18
        FourierParameters -> {1, -1}
19
    ] * FourierTransform[
20
        HeavisidePi[2t], t, \[Omega],
21
        FourierParameters -> {1, -1}
22
    ], {\[Omega], -30, 30}, 
23
    GridLines -> Automatic,
24
    PlotRange -> All
25
] (* 门函数傅里叶变换的乘积, 就是三角波信号的傅里叶变换 *)

3.4 信号的采样过程

xxxxxxxxxx
9
1
F[t_] := Which[
2
    t < -5 Pi, F[t + 10 Pi],
3
    -5 Pi <= t <= 5 Pi, FourierTransform[
4
        HeavisideLambda[2 Subscript[x, temp]],
5
        Subscript[x, temp], t
6
    ],
7
    t > 5 Pi, F[t - 10 Pi]
8
] (* 注意不能使用 F[Subscript[x, temp]] (所以我才用了这么个变量名) *)
9
Plot[F[t], {t, -24 Pi, 24 Pi}]

3.5 系统的频率特性

3.5.1 低通滤波器的频率响应

xxxxxxxxxx
3
1
H = 1/(0.08 (I \[Omega])^2 + 0.4 (I \[Omega]) + 1);
2
Plot[Abs[H], {\[Omega], 0, 50}, GridLines -> Automatic]
3
Plot[Arg[H], {\[Omega], 0, 50}, GridLines -> Automatic]

3.5.2 连续 LTI 系统的频率响应

xxxxxxxxxx
5
1
H = (10 I \[Omega] + 5)/(
2
    (I \[Omega])^3 + 10 (I \[Omega])^2 + 8 I \[Omega] + 5
3
);
4
Plot[Abs[H], {\[Omega], -8, 8}, GridLines -> Automatic]
5
Plot[Arg[H], {\[Omega], -8, 8}, GridLines -> Automatic]

第 4 章 连续复频域分析

平时一直用的 mathematica，这里就不再写代码了.

4.0 复频域常用函数

4.1 拉普拉斯变换

4.2 与傅氏变换关系

4.3 拉普拉斯逆变换

4.4 系统函数零极点图

4.5 系统的复频域分析

4.6 系统的频率响应图

第 5 章 离散时域分析

5.1 常用函数与功能

5.1.1 离散时间序列的绘图

xxxxxxxxxx
1
1
DiscretePlot[DiscreteDelta[n - 1], {n, -2, 2}]

5.1.2 离散序列卷积和的求解

xxxxxxxxxx
15
1
x[t_] := Switch[t,
2
  0, 1,
3
  1, 0,
4
  2, 3,
5
  Default, 0
6
];
7

8
y[t_] := Switch[t,
9
  0, 1,
10
  1, -2,
11
  Default, 0
12
];
13

14
result = DiscreteConvolve[x[m], y[m], m, n];
15
Table[result /. n -> i, {i, 0, 3}]

5.1.3 卷积求解过程的图像绘制

xxxxxxxxxx
5
1
f[n_] := UnitBox[n/7];
2
g[n_] := 2 UnitBox[n/3 + 1];
3
GraphicsRow@
4
 Table[DiscretePlot[{f[x], g[x - i]}, {x, -10, 10}, ImageSize -> 110, 
5
   Ticks -> None, PlotRange -> All], {i, -1, 2}]

5.1.4 十分有趣的交互式动画

xxxxxxxxxx
27
1
DynamicModule[{f, g, h},
2
 f[n_] := (3/4)^n UnitStep[n];
3
 g[n_] := UnitBox[n/5];
4
 h[n_] = DiscreteConvolve[f[m], g[m], m, n];
5
 Manipulate[
6
  GraphicsRow[{DiscretePlot[{f[m], g[n - m]}, {m, -5, 12}, Sequence[
7
     PlotRange -> {-1, 1.5}, Ticks -> {Automatic, {0, 1}}, 
8
      AxesLabel -> Automatic, PlotStyle -> PointSize[0.04], 
9
      PlotLabel -> Text[
10
Style["f[m], g[n-m]", Italic]], Epilog -> {
11
Arrowheads[Medium], 
12
Arrow[{{n, -1}, {n, -0.1}}], 
13
Text[
14
Style["n", Italic], {n + 1.5, -0.9}]}]],
15
    DiscretePlot[f[m] g[n - m], {m, -5, 12}, Sequence[
16
     PlotRange -> {-1, 1.5}, Ticks -> {Automatic, {0, 1}}, 
17
      AxesLabel -> Automatic, PlotStyle -> PointSize[0.04], 
18
      PlotLabel -> Text[
19
Style["f[m]\[CenterDot]g[n-m]", Italic]]]],
20
    DiscretePlot[h[m], {m, -5, n}, Sequence[
21
     PlotRange -> {{-5, 12}, {-1, 4}}, Ticks -> {Automatic, {0, 2, 4}},
22
       AxesLabel -> {
23
Text[
24
Style["n", Italic]], None}, PlotStyle -> PointSize[0.04], 
25
      PlotLabel -> Text[
26
Style["(f*g)[n]", Italic]]]]}, ImageSize -> 480], {{n, 1}, -4, 11, 
27
   1}]]

5.1.5 导出卷积和的演示动画

xxxxxxxxxx
41
1
DynamicModule[{f, g, h},
2
    f[n_] := (3/4)^n UnitStep[n];
3
    g[n_] := UnitBox[n/5];
4
    h[n_] = DiscreteConvolve[f[m], g[m], m, n];
5
    
6
    frames = Table[GraphicsRow[{
7
        DiscretePlot[
8
            {f[m], g[n - m]}, {m, -5, 12}, 
9
            PlotRange -> {-1, 1.5},
10
            Ticks -> {Automatic, {0, 1}}, 
11
            AxesLabel -> Automatic,
12
            PlotStyle -> PointSize[0.04], 
13
            PlotLabel -> Text[Style["f[m], g[n-m]", Italic]], 
14
            Epilog -> {
15
                Arrowheads[Medium],
16
                Arrow[{{n, -1}, {n, -0.1}}], 
17
                Text[Style["n", Italic], {n + 1.5, -0.9}]
18
            } (*Epilog 用于绘制箭头*)
19
        ],
20
        DiscretePlot[
21
            f[m] g[n - m], {m, -5, 12},
22
            PlotRange -> {-1, 1.5},
23
            Ticks -> {Automatic, {0, 1}},
24
            AxesLabel -> Automatic, 
25
            PlotStyle -> PointSize[0.04], 
26
            PlotLabel -> Text[Style["f[m]\[CenterDot]g[n-m]", Italic]]
27
        ], 
28
        DiscretePlot[
29
            h[m], {m, -5, n},
30
            PlotRange -> {{-5, 12}, {-1, 4}},
31
            Ticks -> {Automatic, {0, 2, 4}}, 
32
            AxesLabel -> {Text[Style["n", Italic]], None}, 
33
            PlotStyle -> PointSize[0.04], 
34
            PlotLabel -> Text[Style["(f*g)[n]", Italic]]
35
        ]
36
    }], {n, -4, 11}];
37
    Export[
38
        "image/DiscreteConvolve.gif", frames,
39
        "DisplayDurations" -> 0.5
40
    ];
41
];